ST5226 Spatial Statistics : Geostatistics and Spatial Predictions
Yeo Ming Jie, Jonathan (A0164616A)
2022-11-16
Lecture 7: Introduction to Geostatistics
Loading required packages:
library(sp)
library(gstat)
library(rgdal)
library(RColorBrewer)
library(classInt)Meuse Dataset
Recall the Meuse river bank dataset which contained measurements of the metal concentration along the Meuse river in the Netherlands.
data(meuse)
coordinates(meuse) <- ~x+y # converts meuse into SpatialPointsDataFrame
meuse$lzn <- log(meuse$zinc)
spplot(meuse, "lzn", main="Log-Zinc Measurements",
col.regions = heat.colors(5))
Observe that the zinc concentrations is larger close to the river bank, clear spatial trend in a given direction - zinc concentrations decreasing as we move in a direction perpendicular to the river bank.
Exploratory Variogram Analysis
Checking for existence/strength of spatial correlation from scatter plots of pairs \(Z(s_i)\) and \(Z(s_j)\) grouped according to their separation distances \(h = ||s_i-s_j||\), i.e. h-scatter plots. Consider the plots for the following sets of distance bins:
hscat(log(zinc) ~ 1, meuse, (0:9)*100)
Plots indicate the correlation coefficient \(r\), observe correlation decays with increasing distance.
Next, we plot the sample semivariogram estimator for meuse. Observe the following key features, e.g. nugget, sill \(\approx 0.6\), Range \(\approx 800\).
vgm_log_zinc <- variogram(log(zinc) ~ 1, meuse) # Default in the North-South direction
# ~ 1 defines a single constant predictor, leading to a spatially constant mean coefficient
plot(vgm_log_zinc)
We can also investigate anisotropy i.e. the covariance function \(C(\cdot)\) of the process is a function of both distance and direction, by computing the sample variogram in 4 different directions.
# alpha vector of angles in degrees, clockwise from North
plot(variogram(log(zinc) ~ 1, meuse, alpha=c(0, 45, 90, 135)))
Observe that the variogram grows slowly in the NE direction, but grows most rapidly in the SE direction.
Smoky Mountains Dataset (pH measurements)
Next, consider the dataset containing the pH measurements at locations within the Great Smoky Mountains. Observe from the plot that observations with high pH values lay around the boundary of the region.
smoky_data <- read.table("/Users/jyeo_/Desktop/MSc Statistics Coursework/ST5226 Spatial Statistics/Data/smoky/smoky.dat", header=TRUE, skip=1)
coordinates(smoky_data) <- ~easting + northing # initialize Spatial Points Data Frame
smok <- readOGR("smoky", "smokypoly") # Spatial Polygons## OGR data source with driver: ESRI Shapefile
## Source: "/Users/jyeo_/Desktop/MSc Statistics Coursework/ST5226 Spatial Statistics/Data/smoky", layer: "smokypoly"
## with 1 features
## It has 1 fieldspal <- brewer.pal(4, "RdYlGn")
q4 <- classIntervals(smoky_data$ph, n=4, style="fisher")
spplot(smoky_data, "ph", cuts=q4$brks , col.regions=pal, colorkey=TRUE,
sp.layout=list("sp.polygons", smok), scales=list(draw=TRUE))
Note that the variogram object contains information on the exact values of \(h\) used, along with the nuymber of points used to compute the estimate. Plotting the sample variogram estimator, we obtain the following:
vgm_ph <- variogram(ph ~ 1, smoky_data, cutoff=120, width=10)
head(vgm_ph, n=2)
## np dist gamma dir.hor dir.ver id
## 1 54 6.453213 0.06594352 0 0 var1
## 2 213 15.010015 0.12051268 0 0 var1
plot(vgm_ph)
Available Semivariogram Types
show.vgms(par.strip.text=list(cex=0.6))
In order to get a full list of variogram names used in
gstat, we can type vgm(). The
vgm() function is used to generate variogram models (or add
to existing models).
head(vgm())## short long
## 1 Nug Nug (nugget)
## 2 Exp Exp (exponential)
## 3 Sph Sph (spherical)
## 4 Gau Gau (gaussian)
## 5 Exc Exclass (Exponential class/stable)
## 6 Mat Mat (Matern)Parametric Semivariograms
Suppose we wish to create a Spherical variogram with parameters \(c_0=1\) (nugget), \(c_s=2\) (partial sill) and \(a_s = 10\) (range).
sph1 <- vgm(psill=2, model='Sph', range=10, nugget=1)
class(sph1); sph1## [1] "variogramModel" "data.frame"## model psill range
## 1 Nug 1 0
## 2 Sph 2 10# Plot of a particular semivariogram (Spherical)
show.vgms(sill=2, range=10, models='Sph', nugget=1, max=13)
Alternatively, consider the Gaussian semivariogram model with effective range \(3a_e = 10\)
gauss1 <- vgm(psill=2, model='Gau', range=10, nugget=1)
gauss1## model psill range
## 1 Nug 1 0
## 2 Gau 2 10Adding Semivariograms
The sum of two valid semivariograms is also a valid semivariogram. Consider the addition of the Spherical and Gaussian semivariograms.
gauss2 <- vgm(psill=2, model='Gau', range=10, nugget=0, add.to=sph1)
out <- variogramLine(gauss2, maxdist=15)
plot(out, type='l', main="Sum of Two Variograms")
Fitting Semivariograms using Likelihood-based Methods
fit_v <- fit.variogram(vgm_log_zinc, vgm(1, "Sph", 800, 1))
fit_v## model psill range
## 1 Nug 0.05065923 0.0000
## 2 Sph 0.59060463 896.9976attr(fit_v, "SSErr")## [1] 9.011194e-06plot(vgm_log_zinc, model=fit_v)
Lecture 8: Spatial Prediction (Kriging)
data(meuse.grid)
coordinates(meuse.grid) <- ~x+y
meuse.grid <- as(meuse.grid, "SpatialPixelsDataFrame")
idw.out <- gstat::idw(log(zinc) ~ 1, meuse, meuse.grid, idp=2.5)## [inverse distance weighted interpolation]plot(idw.out, col=heat.colors(100))
head(as.data.frame(idw.out), n=2)## x y var1.pred var1.var
## 1 181180 333740 6.390860 NA
## 2 181140 333700 6.561591 NAspplot(meuse, "lzn", main="Log-Zinc Measurements",
col.regions = heat.colors(4),key.space="right")
vgm_log_zinc <- variogram(log(zinc) ~ 1, meuse)
fit_v <- fit.variogram(vgm_log_zinc, vgm(1, "Sph", 800, 1))
lz.ok <- krige(log(zinc)~1, meuse, meuse.grid, fit_v)## [using ordinary kriging]plot(lz.ok, col=heat.colors(100))
plot(lz.ok['var1.var'])
plot(meuse, add=TRUE, col='grey50')
err.fit <- fit.variogram(vgm_log_zinc, vgm(1, "Sph", 800, Err=1))
krige(log(zinc) ~ 1, meuse, meuse[1,], err.fit)## [using ordinary kriging]## coordinates var1.pred var1.var
## 1 (181072, 333611) 6.884405 0.03648707krige(log(zinc) ~ 1, meuse, meuse[1,], fit_v)## [using ordinary kriging]## coordinates var1.pred var1.var
## 1 (181072, 333611) 6.929517 4.907917e-34vt.fit2 <- fit.variogram.gls(log(zinc) ~ sqrt(dist), data=meuse[1:40,],
model= vgm(1, "Exp", 300, 1) ,trace=FALSE)
lz.uk.gls <- krige(log(zinc)~sqrt(dist), meuse, meuse.grid, vt.fit2)## [using universal kriging]plot(lz.uk.gls, col=heat.colors(100))
Lecture 9: Spatial Random Effects Model (Meuse data)
Required Packages:
library("sp")
library("FRK")
library("gstat")
library("ggplot2")
library("gridExtra")
library("grid")
library("INLA")Construcing Basic Areal Units (BAUs)
data(meuse)
coordinates(meuse) = ~x+y # changes meuse into SpatialPointsDataFrame
# BAUs from hexagonal grid
HexPols_df <- auto_BAUs(manifold = plane(),
cellsize = 100,
type = "hex",
data = meuse,
nonconvex_hull = FALSE)
# BAUs cover the convex hull of observation locations
plot(HexPols_df)
Constructing the basis functions
# Constructing basis functions
G <- auto_basis(manifold = plane(),data=meuse, # manifold - nature of spatial domain, plane() or sphere()
nres = 2,regular=2,prune=0.1,type = "bisquare") ## Warning in auto_basis(manifold = plane(), data = meuse, nres = 2, regular = 2, :
## it's not considered best practice to use prune anymore and is being deprecated.
## Please set to 0 as it will be removed in future versions# nres - number of resolutions (in this case, small and large scale, 2 resolutions)
show_basis(G,ggplot()) + geom_point(data=data.frame(meuse),aes(x,y))## Note: show_basis assumes spherical distance functions when plotting
Fitting the Spatial Random Effects (SRE) Model
set.seed(1)
# formula for SRE model (just an intercept here)
f <- log(zinc) ~ 1
# Call FRK()
S <- FRK(f = f, # formula
data = list(meuse), # list of data objects (just meuse in this case)
basis=G,
BAUs=HexPols_df) ## Assuming fine-scale variation is homoscedastic
## Modelling using 264 basis functions.
## Constructing SRE model...
## Fitting SRE model...
##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|=== | 4%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 24%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%Minimum tolerance reachedsummary(S) # Print out a summary of the returned (fitted) SRE model## $summ_Kdiag
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.4083 0.4083 0.4083 0.5604 0.4083 1.4375
##
## $summ_mueta
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.41307 -0.17122 0.01199 0.02334 0.22023 1.24154
##
## $summ_vareta
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.2429 0.2846 0.3168 0.4089 0.3913 1.4082
##
## $summ_muxi
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
##
## $finescale_var
## [1] "0.0473070019826517"
##
## $reg_coeff
## [1] "5.91504129250237"
##
## attr(,"class")
## [1] "summary.SRE"Prediction
# Get predictions (for the BAUs by default)
GridBAUsPred<- predict(S, obs_fs = FALSE)
# Convert the SpatialPolygonsDataFrame "GridBAUs1" to an ordinary data frame
# For plotting:
Pred_df <- SpatialPolygonsDataFrame_to_df(sp_polys = GridBAUsPred,vars=c("mu","sd"))
# Set coordinate reference system for predictions
# 4326 is the EPSG code for WGS84
# Create plot for prediction means, mu
g1<-ggplot() +
geom_polygon(
aes(x=x, y=y, fill=mu, group=id), # fill color by mu
data = Pred_df)
# Change x and y axis labels
g1<-g1 + xlab("Easting") + ylab("Northing")
# Plot for standard deviation
g2<-ggplot() +
geom_polygon(
aes(x=x, y=y, fill=sd, group=id),
data = Pred_df)
g2<-g2 + xlab("Easting") + ylab("Northing")
# Combine the above two plots on one graph
grid.arrange(g1,g2,nrow=1)
Lecture 10: Spatial Prediction in Model-Based Geostatistics
library("sp")
library("FRK")
library("gstat")
library("ggplot2")
library("gridExtra")
library("grid")
library("INLA")
library("tidyr")##
## Attaching package: 'tidyr'## The following objects are masked from 'package:Matrix':
##
## expand, pack, unpacklibrary("reshape2")##
## Attaching package: 'reshape2'## The following object is masked from 'package:tidyr':
##
## smithslibrary("gtable")
library("dplyr")##
## Attaching package: 'dplyr'## The following object is masked from 'package:gridExtra':
##
## combine## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionPrediction with Covariates
# Predicting at the BAUs when the large scale trend has covariates
# load meuse.grid and convert to SpatialPixelsDataFrame
data("meuse", package="sp")
coordinates(meuse) = ~ x + y
data("meuse.grid", package = "sp")
coordinates(meuse.grid) = ~ x + y
gridded(meuse.grid) = TRUE
# Remove any common fields in the data and the BAUs (recall that in FRK all covariates
# need to be with the BAUs)
meuse$soil <- meuse$dist <- meuse$ffreq <- NULL
# formula for SRE
f <- log(zinc) ~ 1 + sqrt(dist)
# Run FRK
S <- FRK(f = f, # formula
data = list(meuse), # data (just meuse)
BAUs = meuse.grid, # the BAUs
regular = 0) # irregularly placed basis functions## Assuming fine-scale variation is homoscedastic
## Modelling using 374 basis functions.
## Constructing SRE model...
## Fitting SRE model...
##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|=== | 4%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 24%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%Minimum tolerance reached# Predict at all the BAUs
Pred <- predict(S, obs_fs = FALSE)
Pred_df<-SpatialPolygonsDataFrame_to_df(sp_polys = Pred,vars=c("mu","sd"))g1<-ggplot() +
geom_polygon(aes(x=x, y=y, fill=mu,group=id),
data=Pred_df)
g1<-g1 + xlab("Easting") + ylab("Northing")
g2<-ggplot() +
geom_polygon(aes(x=x, y=y, fill=sd, group=id),
data=Pred_df)
g2<-g2 + xlab("Easting") + ylab("Northing")
grid.arrange(g1,g2,nrow=1)
Data and Prediction Locations in Same Data Frame (missing response data)
# Example where both observed data and prediction locations are
# supplied in one data frame, with response at prediction locations
# missing (NA)
# Reload meuse data and assume first 10 observations are missing
data("meuse", package = "sp")
meuse[1:10, "zinc"] <- NA
# Create a complete-data data frame
meuse2 <- subset(meuse, !is.na(zinc))
meuse2 <- meuse2[, c("x", "y", "zinc")]
coordinates(meuse2) <- ~x+y
# Create BAUs around all prediction and observation
# points (after removing data from field)
meuse$zinc <- NULL
coordinates(meuse) <- c("x", "y")
meuse.grid2 <- BAUs_from_points(meuse)
# Run the function FRK with the data frame and created BAUs.
f <- log(zinc) ~ 1 + sqrt(dist)
S <- FRK(f = f,
data = list(meuse2),
BAUs = meuse.grid2,
regular = 0)## Assuming fine-scale variation is homoscedastic
## Modelling using 403 basis functions.
## Constructing SRE model...
## Averaging over polygons...
## Fitting SRE model...
##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|=== | 4%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 24%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|================================================ | 68%
|
|================================================ | 69%Minimum tolerance reachedPred <- predict(S, obs_fs = FALSE)
data(meuse)
plot(Pred$mu[1:10],log(meuse$zinc[1:10]),xlab="Predicted",ylab="Actual",cex.lab=1.5)
abline(c(0,1))
Predictions of Spatial Averages
# Construct model for meuse data with square grid, and consider
# prediction grid with different support
# Initialise seed to ensure reproducibility
set.seed(1)
# Load meuse data and convert to sp object
data("meuse", package = "sp")
coordinates(meuse) = ~ x + y
# Construct the BAUs
GridBAUs1 <- auto_BAUs(manifold = plane(), # we are on the plane
type = "grid", # gridded BAUs (not hex)
cellsize = c(100, 100), # 100m x 100m
data = meuse, # data for boundary
nonconvex_hull = TRUE, # nonconvex boundary
convex = -0.05) # convex parameter (see
# INLA::inla.nonconvex.hull)G <- auto_basis(manifold = plane(), # we are on the plane
data = meuse, # data around which to construct basis
regular = 0, # irregularly placed basis
nres = 3, # three resolutions
type = "bisquare") # bisquare functions
## Plot basis functions
dev.new()
print(show_basis(G) + # main call
xlab("Easting (m)") + # x-label
ylab("Northing (m)") + # y-label
coord_fixed()) # fixed asp. ratio## Note: show_basis assumes spherical distance functions when plotting# formula (just intercept model)
f <- log(zinc) ~ 1
S <- FRK(f = f,
data = list(meuse),
BAUs = GridBAUs1,
basis = G,
average_in_BAU = FALSE) ## Assuming fine-scale variation is homoscedastic
## Modelling using 374 basis functions.
## Constructing SRE model...
## Fitting SRE model...
##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|=== | 4%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 24%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%Minimum tolerance reached# Now look at prediction of averages on a 600m x 600m grid
Pred_regions <- auto_BAUs(manifold = plane(), # we are on the plane
cellsize = c(600, 600), # 600m x 600m
type = "grid", # grid
data = meuse, # meuse data for boundary
convex = -0.05) # convex parameter (see
# INLA::inla.nonconvex.hull)
## Now predict over these regions
Pred <- predict(S, newdata = Pred_regions, # prediction regions
obs_fs = FALSE) # Case 2 in paper (default)
## Plot results
## First converg SpatialPolygonsDataFrame to data frame
Pred_regions_df <- SpatialPolygonsDataFrame_to_df(sp_polys = Pred, # the sp obj
vars = c("mu", "sd")) # the variables to put
# into the data frameg1<-ggplot() +
geom_polygon(aes(x=x, y=y, fill=mu,group=id),
data=Pred_regions_df)
g1<-g1 + xlab("Easting") + ylab("Northing")
g2<-ggplot() +
geom_polygon(aes(x=x, y=y, fill=sd, group=id),
data=Pred_regions_df)
g2<-g2 + xlab("Easting") + ylab("Northing")
grid.arrange(g1,g2,nrow=1)
Data which are Spatial Averages (Areal Data)
# Now consider an example where the observed data do not have
# point support
set.seed(1) # Fix seed
data("meuse.grid", package = "sp") # load meuse data grid
data("meuse", package = "sp") # load meuse data
## Generate observations with large spatial support
meuse_pols <- NULL # initialise
offset <- 150 # half-size of box (150 m)
## for each meuse data point
for (i in 1:nrow(meuse)) {
this_meuse <- meuse[i, ]
meuse_pols <- rbind(meuse_pols, # create a box aroun it of sides 300 m in length
data.frame(x = c(this_meuse$x - offset,
this_meuse$x + offset,
this_meuse$x + offset,
this_meuse$x - offset),
y = c(this_meuse$y - offset,
this_meuse$y - offset,
this_meuse$y + offset,
this_meuse$y + offset),
id = i,
zinc = this_meuse$zinc))
}
## Convert to SpatialPolygonsDataFrame
meuse_pols <- df_to_SpatialPolygons(meuse_pols,coords = c("x", "y"), keys = "id", proj = CRS())
meuse_pols <- SpatialPolygonsDataFrame(meuse_pols,
data.frame(row.names = row.names(meuse_pols),
zinc = meuse$zinc))
coordnames(meuse_pols) <- c("x", "y")
## Convert meuse and meuse.grid to sp object
coordinates(meuse) <- ~x + y
coordinates(meuse.grid) <- ~x+y# Conditionally simulate field on the BAUs and aggregate to polygons
# to generate realistic observations
# First construct the BAUs
GridBAUs <- auto_BAUs(manifold = plane(), # 2D plane
cellsize = c(50, 50), # BAU cellsize
type = "grid", # grid (not hex)
data = meuse, # data around which to create BAUs
convex = -0.05) # border buffer factor
GridBAUs$fs <- 1
## Now fit variogram to meuse data and conditionally simulate
f <- log(zinc) ~ 1
lzn.vgm <- variogram(f, meuse)
lzn.fit <- fit.variogram(lzn.vgm, model = vgm(1, "Sph", 900, 1))
lzn.sim <- krige(formula = f, meuse,
GridBAUs,
model = lzn.fit,
nsim=1, nmax = 30)## drawing 1 GLS realisation of beta...
## [using conditional Gaussian simulation]# Observe field using large spatial support. In the below we run FRK to just generate
# the mapping matrix, which we then use to assign values to the meuse polygons
S <- FRK(f = f, ## Just get out C Matrix
data = list(meuse_pols),
BAUs = GridBAUs,
regular = 0,
n_EM = 1, print_lik = FALSE)## Modelling using 374 basis functions.
## Constructing SRE model...
## SpatialPolygons were provided for the data support. No 'wts' field was found in the BAUs, so all BAUs are assumed to be of equal weight; if this is not the case, set the 'wts' field in the BAUs accordingly.
## Fitting SRE model...
##
|
| | 0%
|
|======================================================================| 100%
## Maximum EM iterations reachedmeuse_pols$zinc <- as.numeric(S@Cmat %*% exp(as.numeric(lzn.sim$sim1)) +
0.1*rnorm(length(meuse_pols)))## Generate the basis functions
G <- auto_basis(manifold = plane(), # 2D plane
data = meuse, # meuse data
nres = 3, # number of resolutions
type = "bisquare", # type of basis function
regular = 0) # place irregularly in domain
## Construct and fit the SRE model
S <- FRK(f = f,
data=list(meuse_pols),
BAUs=GridBAUs,
basis=G)## Modelling using 374 basis functions.
## Constructing SRE model...
## SpatialPolygons were provided for the data support. No 'wts' field was found in the BAUs, so all BAUs are assumed to be of equal weight; if this is not the case, set the 'wts' field in the BAUs accordingly.
## Fitting SRE model...
##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|=== | 4%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 24%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%Minimum tolerance reachedGridBAUs2 <- predict(S, obs_fs = FALSE)## Create data frames for plotting
BAUs_df <- data.frame(GridBAUs2)
Obs_df <- SpatialPolygonsDataFrame_to_df(sp_polys = meuse_pols, # BAUs to convert
vars = c("zinc")) # fields to extract
## Now plot the footprints
g1 <- ggplot() +
geom_path(data = Obs_df,
aes(x, y, group = id),
colour = "black") +
coord_fixed() +
xlab("Easting (m)") + ylab("Northing (m)")
## Now plot the prediction
g2 <- ggplot() + # Use a plain theme
geom_tile(data = BAUs_df, # Draw BAUs
aes(x, y, fill = mu), # Colour <-> Mean
colour = "light grey") + # Border is light grey
scale_fill_distiller(palette = "Spectral", name = "pred") + # Spectral palette
coord_fixed() + # fix aspect ratio
xlab("Easting (m)") + ylab("Northing (m)") # axes labels
## Now plot the prediction error
g3 <- ggplot() + # Similar to above but with s.e.
geom_tile(data = BAUs_df,
aes(x, y, fill = sqrt(var)),
colour = "light grey") +
scale_fill_distiller(palette = "BrBG",
name = "s.e.") +
coord_fixed() +
xlab("Easting (m)") + ylab("Northing (m)")
dev.new()
plot(g1)#library(gtable)
#p2 <- ggplotGrob(g2)
#p3 <- ggplotGrob(g3)
#g <- cbind(p2, p3, size="first")
#g$heights <- unit.pmax(p2$heights, p3$heights)
#grid.newpage()
#grid.draw(g)
grid.arrange(g2, g3, nrow = 1) # Plot
x.seq <- y.seq <- seq(0, 1, length.out = 100)
BAUs <- expand.grid(x = x.seq, y = y.seq)
coordinates(BAUs) = ~ x + y
gridded(BAUs) <- TRUE
BAUs_df <- coordinates(BAUs) %>% as.data.frameset.seed(2020)
# Define some complicated smooth trigonometric field
#f <- function(x, y, a = 0, b = 0, l = 1) exp(-l * sqrt((x - a)^2 + (y - b)^2))
logistic <- function(x, L = 1, k=1, x0 = 0) L / (1 + exp(-k * (x - x0)))
smooth_Y_process <- function(x, y) {
a <- 4 + 2 * sin(5 * x) + 2 * cos(4 * y) + 2 * sin(3 * x * y)+
sin(7 * x) + cos(9 * y) +
sin(12 * x) + cos(17 * y) + sin(14 * x) * cos(16 * y) +
sin(23 * x) + cos(22 * y) + sin(24 * x) * cos(26 * y) +
x + y + x^2 + y^2 + x * y +
sin(10 * pi * x * y) + cos(20 * x * y - 3) + sin(40 * x * y)
logistic(a, x0 = 2, L = 6, k = 0.35)
}
BAUs_df <- BAUs_df %>%
dplyr::mutate(Y = smooth_Y_process(x, y) + rnorm(length(BAUs), mean = 0, sd = 0.2))
## Compute the true mean process, mu, over the BAUs
BAUs_df <- dplyr::mutate(BAUs_df, mu = exp(Y))
## Subsample n of the BAUs to act as observation locations, and simulate data
n <- 750
Poisson_simulated <- sample_n(BAUs_df, n) %>%
dplyr::mutate(Z = rpois(n, lambda = mu)) %>%
dplyr::select(x, y, Z)
## scalar matrix for fine scale variation, and converts to SpatialPixelsDF
BAUs$fs <- rep(1, length(BAUs))
## Convert Poisson_simulated to Spatial* object
coordinates(Poisson_simulated) <- ~ x + yS <- FRK(f = Z ~ 1, data = list(Poisson_simulated),
nres = 3, BAUs = BAUs,
response = "poisson",
link = "log",
K_type = "precision", method = "TMB", est_error = FALSE) ## Modelling using 819 basis functions.
## Constructing SRE model...
## Fitting SRE model...pred_list<- predict(S, type = c("link", "mean"))
## Predictions, uncertainty, and data
plot_list <- plot(S, pred_list$newdata,
labels_from_coordnames = FALSE)g1<-ggplot(BAUs_df) + geom_tile(aes(x, y, fill = mu)) +
scale_fill_distiller(palette = "Spectral") +
labs(x = expression(s[1]), y = expression(s[2])) +
theme_bw() + coord_fixed() +
scale_x_continuous(breaks=c(0.25, 0.75), expand = c(0, 0)) +
scale_y_continuous(breaks=c(0.25, 0.75), expand = c(0, 0)) +
theme(axis.text = element_text(size = 16),
axis.title = element_text(size = 19),
legend.text = element_text(size = 16),
plot.title = element_text(hjust = 0.5, size = 19))
g1
g2<-plot_list$Z + labs(title = plot_list$Z$labels$fill) +
labs(fill = "", colour = "") +
scale_x_continuous(breaks=c(0.25, 0.75), expand = c(0, 0)) +
scale_y_continuous(breaks=c(0.25, 0.75), expand = c(0, 0)) +
theme(axis.text = element_text(size = 16),
axis.title = element_text(size = 19),
legend.text = element_text(size = 16),
plot.title = element_text(hjust = 0.5, size = 19))
g2
g3<-plot_list$p_Y + labs(title = plot_list$p_Y$labels$fill) +
labs(fill = "", colour = "") +
scale_x_continuous(breaks=c(0.25, 0.75), expand = c(0, 0)) +
scale_y_continuous(breaks=c(0.25, 0.75), expand = c(0, 0)) +
theme(axis.text = element_text(size = 16),
axis.title = element_text(size = 19),
legend.text = element_text(size = 16),
plot.title = element_text(hjust = 0.5, size = 19))
g3